#include "decoding.h"
#include <iostream>
#include <vector>
#include <thread>
#include <chrono>
#include <mutex>
#include <unistd.h>
#include <sys/sysinfo.h>

using namespace std;
// Вектора дублируются при возвращении из функции (доп нагрузка)

Approx::Approx()
{
    this->poly_data = nullptr;
    this->lineal_data = nullptr;
}

Lineal_data::Lineal_data()
{
    this->pair_count = 0;
    this->borders = nullptr;
    this->koefs = nullptr;
    this->ad_data = nullptr;
}

double calc_polynom_point_vector(vector<double> *polynom_values, double x)
{
    double res = 0;
    for(int i = 0; i < (*polynom_values).size(); i++)
    {
        res += (*polynom_values)[i] * pow(x, i);
    }
    return res;
}

template <typename Ty>
double Approx::calc_polynom_point(Ty x)
{
    if(this->poly_data == nullptr)
    {
        throw exception();
    }
    double res = 0;
    for(int i = 0; i < this->poly_data->size; i++)
    {
        res += this->poly_data->poly[i] * pow(x, i);
    }
    return res;
}

void Approx::calc_polynom_func(double *x_values, double *output, int x_size)
{
    for(int i = 0; i < x_size; i++)
    {
        output[i] = calc_polynom_point(x_values[i]);
    }
}

void Approx::get_poly(double *output)
{
    if(this->poly_data == nullptr)
    {
        return;
    }
    this->poly_data->copy_poly(output);
}

int Approx::get_poly_size()
{
    return this->poly_data->size;
}

double Approx::get_mistake()
{
    return this->poly_data->mistake;
}

void Approx::set_poly(double *poly, int size, double mistake)
{
    this->poly_data = new Poly_data;
    this->poly_data->size = size;
    this->poly_data->mistake = mistake;
    this->poly_data->poly = new double[size];
    for(int i = 0; i < size; i++)
    {
        this->poly_data->poly[i] = poly[i];
    }
}

void Approx::get_koefs(double *output)
{
    for(int i = 0; i < this->lineal_data->pair_count; i++)
    {
        output[i * 2] = this->lineal_data->koefs[i * 2];
        output[i * 2 + 1] = this->lineal_data->koefs[i * 2 + 1];
    }
}

void Approx::set_koefs(double *koefs, double *borders, int borders_count)
{
    this->lineal_data = new Lineal_data;
    this->lineal_data->borders = new double[borders_count];
    this->lineal_data->koefs = new double[borders_count * 2];
    this->lineal_data->pair_count = borders_count;
    for(int i = 0; i < borders_count; i++)
    {
        this->lineal_data->koefs[i * 2] = koefs[i * 2];
        this->lineal_data->koefs[i * 2 + 1] = koefs[i * 2 + 1];
        this->lineal_data->borders[i] = borders[i];
    }
}

void Approx::lineal_pices(double *x_values, double *y_values, int input_size)
{
    if(x_values == NULL || y_values == NULL)
    {
        this->lineal_data = nullptr;
        return;
    }
    this->lineal_data = new Lineal_data;
    this->lineal_data->copy_to_borders(x_values, input_size);
    this->lineal_data->koefs = new double[input_size * 2];
    this->lineal_data->pair_count = input_size;
    for(int i = 0; i < input_size - 1; i++)
    {
        this->lineal_data->koefs[i*2+1] = (y_values[i+1] - y_values[i]) / (x_values[i+1] - x_values[i]);
        this->lineal_data->koefs[i*2] = y_values[i] - this->lineal_data->koefs[i*2+1] * x_values[i];
    }
    this->lineal_data->koefs[input_size * 2 - 1] = this->lineal_data->koefs[input_size * 2 - 3];
    this->lineal_data->koefs[input_size * 2 - 2] = this->lineal_data->koefs[input_size * 2 - 4];
}

int Poly_data::copy_poly(double *output)
{
    for(int i = 0; i < this->size; i++)
    {
        output[i] = this->poly[i];
    }
    return this->size;
}

void Lineal_data::slise_size()
{
    this->ad_data->slise_size = this->borders[1] - this->borders[0];
    for(int i = 2; i < this->pair_count; i++)
    {
        if(this->ad_data->slise_size > this->borders[i] - this->borders[i - 1])
        {
            this->ad_data->slise_size = this->borders[i] - this->borders[i - 1];
        }
    }
}

void Lineal_data::slise_count()
{
    this->ad_data->slise_count = (this->borders[this->pair_count - 1] - this->borders[0]) / this->ad_data->slise_size;
}

void Lineal_data::lineal_split()
{
    this->ad_data->slises = new int[this->ad_data->slise_count];
    int borders_ind = 0;
    for(int i = 0; i < this->ad_data->slise_count; i++)
    {
        double pos = this->ad_data->slise_size * i + this->borders[0];
        if(borders_ind < this->pair_count + 1 and pos > this->borders[borders_ind + 1])
        {
            borders_ind++;
        }
        this->ad_data->slises[i] = borders_ind;
    }
}

void Lineal_data::fill_ad_data()
{
    this->ad_data = new Additional_data;
    this->slise_size();
    this->slise_count();
    this->lineal_split();
}

void Lineal_data::copy_to_borders(double *values, int count)
{
    this->borders = new double[count];
    for(int i = 0; i < count; i++)
    {
        this->borders[i] = values[i];
    }
}

template <typename Ty>
double Approx::calc_split_lineal_point(Ty point)
{
    if(this->lineal_data == nullptr)
    {
        throw exception();
    }
    if(this->lineal_data->ad_data == nullptr)
    {
        this->lineal_data->fill_ad_data();
    }
    int koef_ind = 0;
    if(point > this->lineal_data->borders[0])
    {
        int split_ind = (point - this->lineal_data->borders[0]) / this->lineal_data->ad_data->slise_size;
        if(split_ind >= this->lineal_data->ad_data->slise_count)
        {
            koef_ind = this->lineal_data->pair_count - 1;
        }
        else if(point > this->lineal_data->borders[this->lineal_data->ad_data->slises[split_ind] + 1])
        {
            koef_ind = this->lineal_data->ad_data->slises[split_ind] + 1;
        }
        else
        {
            koef_ind = this->lineal_data->ad_data->slises[split_ind];
        }
    }
    return this->lineal_data->koefs[koef_ind * 2] + this->lineal_data->koefs[koef_ind * 2 +1] * point;
}

void Approx::calc_split_lineal(double *x_values, double *output, int x_size)
{
    for(int i = 0; i < x_size; i++)
    {
        output[i] = calc_split_lineal_point(x_values[i]);
    }
}
int binary_step(double x, double *borders, int left, int right)
{
    int mid = (left + right) / 2;
    if(x >= borders[mid])
    {
        if(mid == right or x <= borders[mid + 1])
            return mid;
        return binary_step(x, borders, mid + 1, right);
    }
    else
    {
        if(mid == left)
            return mid;
        if(x >= borders[mid - 1])
            return mid - 1;
        return binary_step(x, borders, left, mid - 1);
    }
}

int binary_search(double x, double *borders, int borders_count)
{
    return binary_step(x, borders, 0, borders_count - 1);
}

template <typename Ty>
double Approx::calc_lineal_point(Ty point)
{
    if(this->lineal_data == nullptr)
    {
        throw exception();
    }
    int koef_ind = binary_search(point, this->lineal_data->borders, this->lineal_data->pair_count);
    return this->lineal_data->koefs[koef_ind * 2] + this->lineal_data->koefs[koef_ind * 2 +1] * point;
}

void Approx::calc_lineal_pices(double *x_values, double *output, int x_size)
{
    for(int i = 0; i < x_size; i++)
    {
        output[i] = calc_lineal_point(x_values[i]);
    }
}

double vector_mid(vector<double> *v)
{
    double mid = 0;
    for (int i = 0; i < (*v).size(); i++)
    {
        mid += (*v)[i];
    }
    return mid / (*v).size();
}

double difference(vector<double> *x_values, vector<double> *y_values, vector<double> *polynom_values)
{
    double res = 0;
    for(int i = 0; i < (*x_values).size(); i++)
    {
        double y = 0;
        for(int j = 0; j < (*polynom_values).size(); j++)
        {
            y += (*polynom_values)[j] * pow((*x_values)[i], j);
        }
        res += abs(y - (*y_values)[i]);
    }
    return res;
}

void Approx::approx_mistake(double *x_values, double *y_values, int gh_size)
{
    if (this->poly_data == nullptr)
    {
        return;
    }
    double res = 0;
    double *y2_values = new double[gh_size];
    calc_polynom_func(x_values, y2_values, gh_size);
    double max1 = y_values[0];
    double min1 = y_values[0];
    for(int i = 0; i < gh_size; i++)
    {
        if(y_values[i] > max1)
            max1 = y_values[i];
        else if (y_values[i] < min1)
            min1 = y_values[i];
    }
    if (max1 == min1)
    {
        approx_mistake_bigzero(x_values, y_values, gh_size);
        return;
    }
    for(int i = 0; i < gh_size; i++)
    {
        res += abs(y2_values[i] - y_values[i]);
    }
    delete[] y2_values;
    this->poly_data->mistake = res / gh_size / (max1 - min1) * 100;
}

void Approx::approx_mistake_bigzero(double *x_values, double *y_values, int gh_size)
{
    if (this->poly_data == nullptr)
    {
        return;
    }
    double res = 0;
    double *y2_values = new double[gh_size];
    calc_polynom_func(x_values, y2_values, gh_size);
    int n = gh_size;
    for(int i = 0; i < gh_size; i++)
    {
        if(y_values[i] == 0)
            n -= 1;
        else
            res += abs((y2_values[i] - y_values[i]) / y_values[i]);
    }
    delete[] y2_values;
    if(n)
    {
        this->poly_data->mistake = res / n * 100;
        return;
    }
    for(int i = 0; i < gh_size; i++)
    {
        if(y2_values[i] != 0)
        {
            this->poly_data->mistake = 100;
            return;
        }
    }
    this->poly_data->mistake = 0;
}

vector<vector<double>> fill_matrix(vector<double> *x_values, int max_pow)
{
    int size = (*x_values).size();
    vector<vector<double>> matrix(max_pow);
    for(int i = 0; i < max_pow; i++)
    {
        matrix[i].resize(max_pow);
        for(int j = 0; j < max_pow; j++)
        {
            matrix[i][j] = 0;
            for(int k = 0; k < size; k++)
            {
                matrix[i][j] += pow((*x_values)[k], i + j);
            }
        }
    }
    return matrix;
}

void fill_b_vector(vector<double> *vector_b, vector<double> *x_values, vector<double> *y_values, int max_pow)
{
    for(int i = 0; i < max_pow; i++)
    {
        for(int k = 0; k < (*x_values).size(); k++)
        {
            (*vector_b)[i] += pow((*x_values)[k], i) * (*y_values)[k];
        }
    }
}

struct fill_cells_thread_interface
{
public:
    vector<vector<double>> *matrix;
    vector<double> *x_values;
    vector<double> cell;
    std::mutex have_task;
    std::mutex *need_task;
    int *free_count;
    bool *have_job;
    bool free_flag;
    int data_size;

    void init(vector<vector<double>> *matrix, vector<double> *x_values, std::mutex *need_task, bool *have_job, int *free_count)
    {
        this->matrix = matrix;
        this->x_values = x_values;
        this->need_task = need_task;
        this->have_job = have_job;
        this->free_count = free_count;
        free_flag = true;
        have_task.lock();
        cell.resize(2);
        data_size = (*x_values).size();
    }
};

void fill_cells_thread(fill_cells_thread_interface *interface)
{
    int i = -1;
    int j = -1;
    int counter = 0;
    while(*interface->have_job)
    {
        (*interface->need_task).unlock();
        interface->have_task.lock();
        if(i == interface->cell[0] && j == interface->cell[1])
            break;
        i = interface->cell[0];
        j = interface->cell[1];
        double res = 0;
        for(int k = 0; k < interface->data_size; k++)
        {
            res += pow((*interface->x_values)[k], i + j);
        }
        (*interface->matrix)[i][j] = res;
        counter++;
        (*interface->free_count)++;
        interface->free_flag = true;
    }
}

vector<vector<double>> fill_matrix_par(vector<double> *x_values, int max_pow, int thr_count)
{
    int pos = 0;
    if(thr_count >= max_pow * max_pow)
        thr_count = max_pow * max_pow - 1;
    int free_count = thr_count;
    std::mutex need_task;
    vector<fill_cells_thread_interface> thr_interfaces(thr_count);
    bool have_job = true;
    vector<vector<double>> matrix(max_pow);
    for(int i = 0; i < thr_count; i++)
    {
        thr_interfaces[i].init(&matrix, x_values, &need_task, &have_job, &free_count);
    }
    vector<std::thread*> threads(thr_count);
    for(int i = 0; i < thr_count; i++)
    {
        threads[i] = new std::thread(fill_cells_thread, &thr_interfaces[i]);
    }
    int counter = 0;
    for(int i = 0; i < max_pow; i++)
    {
        matrix[i].resize(max_pow);
        for(int j = 0; j < max_pow; j++)
        {
            matrix[i][j] = 0;
            while(true)
            {
                if(not free_count)
                {
                    need_task.lock();
                }
                if(thr_interfaces[pos].free_flag)
                {
                    counter++;
                    thr_interfaces[pos].free_flag = false;
                    free_count--;
                    thr_interfaces[pos].cell[0] = i;
                    thr_interfaces[pos].cell[1] = j;
                    thr_interfaces[pos].have_task.unlock();
                    pos = (pos + 1) % thr_count;
                    break;
                }
                pos = (pos + 1) % thr_count;
            }
        }
    }
    have_job = false;
    for(int i = 0; i < thr_count; i++)
    {
        thr_interfaces[i].have_task.unlock();
    }
    for(int i = 0; i < thr_count; i++)
    {
        (*threads[i]).join();
    }
    return matrix;
}

vector<vector<double>> fill_matrix(double *x_values, int x_size, int max_pow)
{
    vector<vector<double>> matrix(max_pow);
    for(int i = 0; i < max_pow; i++)
    {
        matrix[i].resize(max_pow);
        for(int j = 0; j < max_pow; j++)
        {
            matrix[i][j] = 0;
            for(int k = 0; k < x_size; k++)
            {
                matrix[i][j] += pow(x_values[k], i + j);
            }
        }
    }
    return matrix;
}

void Approx::least_squares_arr(double *x_values, double *y_values, int input_size, int output_pow)
{
    if(x_values == NULL || y_values == NULL)
    {
        this->poly_data = nullptr;
        return;
    }
    this->poly_data = new Poly_data;
    output_pow++;
    this->poly_data->size = output_pow;
    this->poly_data->poly = new double[output_pow];
    for(int i = 0; i < output_pow; i++)
    {
        this->poly_data->poly[i] = 0;
        for(int k = 0; k < input_size; k++)
        {
            this->poly_data->poly[i] += pow(x_values[k], i) * y_values[k];
        }
    }
    vector<vector<double>> matrix_a;
    matrix_a = fill_matrix(x_values, input_size, output_pow);
    for(int i = 0; i < output_pow; i++)
    {
        int not_zero = -1;
        for(int j = i; j < output_pow; j++)
        {
            if(matrix_a[i][j] != 0)
            {
                not_zero = j;
                break;
            }
        }
        if(not_zero == -1)
            continue;
        if(i != not_zero)
        {
            vector<double> tmp_m = matrix_a[i];
            matrix_a[i] = matrix_a[not_zero];
            matrix_a[not_zero] = tmp_m;

            double tmp_v = this->poly_data->poly[i];
            this->poly_data->poly[i] = this->poly_data->poly[not_zero];
            this->poly_data->poly[not_zero] = tmp_v;
        }
        this->poly_data->poly[i] /= matrix_a[i][i];
        for(int j = i + 1; j < output_pow; j++)
        {
            matrix_a[i][j] /= matrix_a[i][i];
        }
        matrix_a[i][i] = 1;
        for(int j = i + 1; j < output_pow; j++)
        {
            if(matrix_a[j][i] == 0)
                continue;
            this->poly_data->poly[j] -= matrix_a[j][i] * this->poly_data->poly[i];
            for(int k = i + 1; k < output_pow; k++)
            {
                matrix_a[j][k] -= matrix_a[j][i] * matrix_a[i][k];
            }
            matrix_a[j][i] = 0;
        }
    }
    for(int i = output_pow - 1; i >= 0; i--)
    {
        if(matrix_a[i][i] == 0)
        {
            this->poly_data->poly[i] = 0;
            continue;
        }
        for(int j = i - 1; j >= 0; j--)
        {
            this->poly_data->poly[j] -= matrix_a[j][i] * this->poly_data->poly[i];
        }
    }

    approx_mistake(x_values, y_values, input_size);
}

struct pascal
{
    vector<vector<double>> triangle;

    void init(int size)
    {
        triangle.resize(size);
        triangle[0].resize(1);
        triangle[0][0] = 1;

        for(int i = 1; i < size; i++)
        {
            triangle[i].resize(i + 1);
            triangle[i][0] = 1;
            for(int j = 1; j < i; j++)
            {
                triangle[i][j] = triangle[i-1][j-1] + triangle[i-1][j];
            }
            triangle[i][i] = 1;
        }
    }
};

struct resizer
{
    double mult;
    double shift;
    int size;
    vector<double> new_x;
    vector<double> new_y;

    void init_resize(vector<double> *x_values, vector<double> *y_values)
    {
        size = (*x_values).size();
        new_x.resize(size);
        new_y.resize(size);

        mult = 2 / ((*x_values)[size - 1] - (*x_values)[0]);
        shift = (*x_values)[0] * mult + 1;

        for(int i = 0; i < size; i++)
        {
            new_x[i] = (*x_values)[i] * mult - shift;
            new_y[i] = (*y_values)[i] * mult;
        }
    }

    void init_resize(double *x_values, double *y_values, int size)
    {
        this->size = size;
        new_x.resize(size);
        new_y.resize(size);

        mult = 2 / (x_values[size - 1] - x_values[0]);
        shift = x_values[0] * mult + 1;

        for(int i = 0; i < size; i++)
        {
            new_x[i] = x_values[i] * mult - shift;
            new_y[i] = y_values[i] * mult;
        }
    }

    void resize_poly(vector<double> *pol, pascal *pasc)
    {
        for(int i = 1; i < (*pol).size(); i++)
            for(int j = 0; j < i; j++)
            {
                (*pol)[j] += (*pol)[i] * pasc->triangle[i][j] * pow(-shift, i - j);
            }

        for(int i = 0; i < (*pol).size(); i++)
        {
            (*pol)[i] *= (1 / mult) * pow(mult, i);
        }
    }
};

vector<double> poly_mult_poly(vector<double> *poly1, vector<double> *poly2)
{
    int size1 = (*poly1).size();
    int size2 = (*poly2).size();
    vector<double> res(size1 + size2 - 1);
    for(int i = 0; i < res.size(); i++)
    {
        res[i] = 0;
    }
    for(int i = 0; i < size1; i++)
        for(int j = 0; j < size2; j++)
        {
            res[i + j] += (*poly1)[i] * (*poly2)[j];
        }
    return res;
}

double integral(double a, double b, vector<double> *poly)
{
    (*poly).push_back(0);
    for(int i = (*poly).size() - 2; i >= 0; i--)
    {
        (*poly)[i+1] = (*poly)[i] / (i + 1);
    }
    (*poly)[0] = 0;
    double y_a = calc_polynom_point_vector(poly, a);
    double y_b = calc_polynom_point_vector(poly, b);
    return y_a - y_b;
}

vector<double> poly_mult_lin(vector<double> *P_i, double k, double b)
{
    vector<double> poly2 = {b, k};
    return poly_mult_poly(P_i, &poly2);
}

vector<vector<double>> get_legendre_polynoms(int max_poly, pascal *pasc)
{
    vector<vector<double>> res(max_poly + 1);
    res[0] = {1};
    res[1] = {0, 1};
    for(int p = 2; p <= max_poly; p++)
    {
        res[p].resize(p + 1);
        for(int i = 0; i <= p; i++)
        {
            res[p][i] = 0;
        }
        vector<double> before_pr(p + 1);
        for(int i = 0; i <= p; i++)
        {
            before_pr[i] = pasc->triangle[p][i] * pow(-1, p - i);
        }

        int fac = 1;
        for(int i = 2; i <= p; i++)
        {
            fac *= i;
        }
        double mult = 1.0/(fac * pow(2, p));

        for(int i = p; i*2 >= p; i--)
        {
            double cur = before_pr[i];
            int p2 = i*2;
            for(int j = p2; j > p2 - p; j--)
            {
                cur *= j;
            }
            res[p][p2 - p] = cur * mult;
        }
    }

    return res;
}

vector<double> legendre_approx(resizer resiz, int max_pow)
{
    //P0 P1 P2 ...
    vector<double> res(max_pow + 1);
    pascal pasc;
    pasc.init(max_pow + 1);
    for(int i = 0; i <= max_pow; i++)
    {
        res[i] = 0;
    }
    vector<vector<double>> P = get_legendre_polynoms(max_pow, &pasc);
    vector<vector<double>> lineal = approx::lineal_pices(&resiz.new_x, &resiz.new_y);
    for(int i = 0; i <= max_pow; i++)
    {
        double c = 0;
        for(int j = 1; j < resiz.size; j++)
        {
            vector<double> tmp = poly_mult_lin(&P[i], lineal[1][j-1], lineal[0][j-1]);
            c += ((2 * i + 1) * 0.5) * integral(resiz.new_x[j], resiz.new_x[j - 1], &tmp);
        }
        for(int j = 0; j < P[i].size(); j++)
        {
            res[j] += c * P[i][j];
        }
    }
    //return res;
    resiz.resize_poly(&res, &pasc);
    return res;
}

void Approx::legendre_approx_arr(double *x_values, double *y_values, int input_size, int output_pow)
{
    if(x_values == NULL || y_values == NULL)
    {
        this->poly_data = nullptr;
        return;
    }
    this->poly_data = new Poly_data;
    this->poly_data->size = output_pow;
    this->poly_data->poly = new double[output_pow];
    resizer resiz;
    resiz.init_resize(x_values, y_values, input_size);
    vector<double> tmp = legendre_approx(resiz, output_pow);
    for(int i = 0; i <= output_pow; i++)
    {
        this->poly_data->poly[i] = tmp[i];
    }
    approx_mistake(x_values, y_values, input_size);
}

void Approx::legendre_approx_auto_pow(double *x_values, double *y_values, int input_size)
{
    if(x_values == NULL || y_values == NULL)
    {
        this->poly_data = nullptr;
        return;
    }
    legendre_approx_arr(x_values, y_values, input_size, 1);
    for(int next_size = 3; next_size <= 8; next_size++)
    {
        Approx new_approx;
        new_approx.legendre_approx_arr(x_values, y_values, input_size, next_size - 1);
        if(new_approx.poly_data->mistake  / (input_size - new_approx.poly_data->size) < this->poly_data->mistake  / (input_size - this->poly_data->size))
        {
            this->poly_data = new_approx.poly_data;
        }
        else
        {
            return;
        }
    }
}

void Approx::least_squares_auto_pow(double *x_values, double *y_values, int input_size)
{
    if(x_values == NULL || y_values == NULL)
    {
        this->poly_data = nullptr;
        return;
    }
    least_squares_arr(x_values, y_values, input_size, 1);
    for(int next_size = 3; next_size <= 8; next_size++)
    {
        Approx new_approx;
        new_approx.least_squares_arr(x_values, y_values, input_size, next_size - 1);
        if(new_approx.poly_data->mistake  / (input_size - new_approx.poly_data->size) < this->poly_data->mistake  / (input_size - this->poly_data->size))
        {
            this->poly_data = new_approx.poly_data;
        }
        else
        {
            return;
        }
    }
}
template double Approx::calc_polynom_point(char point);
template double Approx::calc_polynom_point(unsigned char point);
template double Approx::calc_polynom_point(short point);
template double Approx::calc_polynom_point(unsigned short point);
template double Approx::calc_polynom_point(int point);
template double Approx::calc_polynom_point(unsigned int point);
template double Approx::calc_polynom_point(long long point);
template double Approx::calc_polynom_point(unsigned long long point);
template double Approx::calc_polynom_point(double point);
template double Approx::calc_polynom_point(float point);

template double Approx::calc_lineal_point(char point);
template double Approx::calc_lineal_point(unsigned char point);
template double Approx::calc_lineal_point(short point);
template double Approx::calc_lineal_point(unsigned short point);
template double Approx::calc_lineal_point(int point);
template double Approx::calc_lineal_point(unsigned int point);
template double Approx::calc_lineal_point(long long point);
template double Approx::calc_lineal_point(unsigned long long point);
template double Approx::calc_lineal_point(double point);
template double Approx::calc_lineal_point(float point);

template double Approx::calc_split_lineal_point(char point);
template double Approx::calc_split_lineal_point(unsigned char point);
template double Approx::calc_split_lineal_point(short point);
template double Approx::calc_split_lineal_point(unsigned short point);
template double Approx::calc_split_lineal_point(int point);
template double Approx::calc_split_lineal_point(unsigned int point);
template double Approx::calc_split_lineal_point(long long point);
template double Approx::calc_split_lineal_point(unsigned long long point);
template double Approx::calc_split_lineal_point(double point);
template double Approx::calc_split_lineal_point(float point);

namespace approx{
// !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
vector<double> least_squares(vector<double> *x_values, vector<double> *y_values, int max_pow, int thr_count)
{
    thr_count = 1;
    if(thr_count < 1)
    {
        thr_count = get_nprocs();
    }
    max_pow++;
    vector<double> vector_b(max_pow);
    thread thr(fill_b_vector, &vector_b, x_values, y_values, max_pow);
    vector<vector<double>> matrix_a;
    if(thr_count > 1)
    {
        matrix_a = fill_matrix_par(x_values, max_pow, thr_count);
    }
    else
    {
        matrix_a = fill_matrix(x_values, max_pow);
    }
    thr.join();
    for(int i = 0; i < max_pow; i++)
    {
        int not_zero = -1;
        for(int j = i; j < max_pow; j++)
        {
            if(matrix_a[i][j] != 0)
            {
                not_zero = j;
                break;
            }
        }
        if(not_zero == -1)
            continue;
        if(i != not_zero)
        {
            vector<double> tmp_m = matrix_a[i];
            matrix_a[i] = matrix_a[not_zero];
            matrix_a[not_zero] = tmp_m;

            double tmp_v = vector_b[i];
            vector_b[i] = vector_b[not_zero];
            vector_b[not_zero] = tmp_v;
        }
        vector_b[i] /= matrix_a[i][i];
        for(int j = i + 1; j < max_pow; j++)
        {
            matrix_a[i][j] /= matrix_a[i][i];
        }
        matrix_a[i][i] = 1;
        for(int j = i + 1; j < max_pow; j++)
        {
            if(matrix_a[j][i] == 0)
                continue;
            vector_b[j] -= matrix_a[j][i] * vector_b[i];
            for(int k = i + 1; k < max_pow; k++)
            {
                matrix_a[j][k] -= matrix_a[j][i] * matrix_a[i][k];
            }
            matrix_a[j][i] = 0;
        }
    }
    for(int i = max_pow - 1; i >= 0; i--)
    {
        if(matrix_a[i][i] == 0)
        {
            vector_b[i] = 0;
            continue;
        }
        for(int j = i - 1; j >= 0; j--)
        {
            vector_b[j] -= matrix_a[j][i] * vector_b[i];
        }
    }
    return vector_b;
}

vector<vector<double>> lineal_pices(vector<double> *x_values, vector<double> *y_values)
{
    vector<vector<double>> res(2);
    int res_size = (*x_values).size();
    res[0].resize(res_size);
    res[1].resize(res_size);

    for(int i = 0; i < res_size - 1; i++)
    {
        res[1][i] = ((*y_values)[i+1] - (*y_values)[i]) / ((*x_values)[i+1] - (*x_values)[i]);
        res[0][i] = (*y_values)[i] - res[1][i] * (*x_values)[i];
    }
    res[1][res_size - 1] = res[1][res_size - 2];
    res[0][res_size - 1] = res[0][res_size - 2];

    return res;
}

vector<double> legendre_approx(vector<double> *x_values, vector<double> *y_values, int max_pow)
{
    resizer resiz;
    resiz.init_resize(x_values, y_values);
    return legendre_approx(resiz, max_pow);
}

double squared_deviation(vector<double> *v)
{
    double mid = vector_mid(v);
    double res = 0;
    for (int i = 0; i < (*v).size(); i++)
    {
        res += pow((*v)[i] - mid, 2);
    }
    return sqrt(res / (*v).size());
}

double correlation(vector<double> *x_values, vector<double> *y_values, vector<double> *polynom_values)
{
    double res = 0;
    double y_mid = vector_mid(y_values);
    vector<double> y2_values = calc_polynom_func(polynom_values, x_values);
    double y2_mid = vector_mid(&y2_values);
    for(int i = 0; i < (*y_values).size(); i++)
    {
        res += ((*y_values)[i] - y_mid) * (y2_values[i] - y2_mid);
    }
    res /= (*y_values).size() * squared_deviation(y_values) * squared_deviation(&y2_values);
    return res;
}

double approx_mistake(vector<double> *x_values, vector<double> *y_values, vector<double> *polynom_values)
{
    double res = 0;
    vector<double> y2_values = calc_polynom_func(polynom_values, x_values);
    double max1 = (*y_values)[0];
    double min1 = (*y_values)[0];
    for(int i = 0; i < (*y_values).size(); i++)
    {
        if((*y_values)[i] > max1)
            max1 = (*y_values)[i];
        else if ((*y_values)[i] < min1)
            min1 = (*y_values)[i];
    }
    for(int i = 0; i < (*y_values).size(); i++)
    {
        res += abs(y2_values[i] - (*y_values)[i]);
    }
    return res / (*y_values).size() / (max1 - min1) * 100;
}

vector<double> calc_polynom_func(vector<double> *polynom_values, vector<double> *x_values)
{
    vector<double> res((*x_values).size());
    for(int i = 0; i < (*x_values).size(); i++)
    {
        res[i] = calc_polynom_point_vector(polynom_values, (*x_values)[i]);
    }
    return res;
}

double calc_polynom_point(double *poly, int size, double point)
{
    double res = 0;
    for(int i = 0; i < size; i++)
    {
        res += poly[i] * pow(point, i);
    }
    return res;
}

double calc_lineal_point(double *koefs, double *borders, int borders_count, double point)
{
    int koef_ind = binary_search(point, borders, borders_count);
    return koefs[koef_ind * 2] + koefs[koef_ind * 2 +1] * point;
}

}

namespace MODGH {

    void x5(double *input_output_x, int x_size, double k1, double k2, double k3)
    {
        if(k1 == 0)
            k1 = 1;
        if(k2 == 0)
            k2 = 1;
        if(k3 == 0)
            k3 = 1;

        for(int i = 0; i < x_size; i++)
        {
            input_output_x[i] *= k1 * k2 * k3;
        }
    }

    void x6(double *input_output_x, int x_size, double c1, double c2, double c3)
    {
        for(int i = 0; i < x_size; i++)
        {
            input_output_x[i] += c1 + c2 + c3;
        }
    }

    void x7(double *input_output_x, int x_size, double g1, double g2, double g3)
    {
        if(g1 == 0)
            g1 = 1;
        if(g2 == 0)
            g2 = 1;
        if(g3 == 0)
            g3 = 1;

        for(int i = 0; i < x_size; i++)
        {
            input_output_x[i] /= g1 * g2 * g3;
        }
    }

    void y1(double *input_output_y, int y_size, double Kp)
    {
        if(Kp == 0)
            Kp = 1;

        for(int i = 0; i < y_size; i++)
        {
            input_output_y[i] /= Kp;
        }
    }

    void y2(double *input_output_y, int y_size, double Dp)
    {
        for(int i = 0; i < y_size; i++)
        {
            input_output_y[i] += Dp;
        }
    }

    void y3(double *input_output_y, int y_size, double nn)
    {
        for(int i = 0; i < y_size; i++)
        {
            input_output_y[i] *= nn;
        }
    }

    void y4(double *input_output_y, int y_size, double d, double k)
    {
        if(k == 0)
            k = 1;
        for(int i = 0; i < y_size; i++)
        {
            input_output_y[i] = (input_output_y[i] - d) / k;
        }
    }

    void x1(double *input_output_x, int x_size, double Ub)
    {
        if(Ub == 0)
            Ub = 1;

        for(int i = 0; i < x_size; i++)
        {
            input_output_x[i] = input_output_x[i] / Ub * 100;
        }
    }

    void y5(double *output_2x, double *output_2y, double U0, double Ub, double k)
    {
        if(k == 0)
            k = 1;
        output_2x[0] = 0;
        output_2x[1] = 1;
        output_2y[0] = -U0 * k;
        output_2y[1] = (Ub - U0) * k;
    }

    void xy5(double *output_30x, double *output_30y, double L, double r, double R)
    {
        if(L == 0)
            L = 1;
        if(r == 0)
            r = 1;
        if(R == 0)
            R = 1;
        output_30y[0] = 0;
        for(int i = 1; i < 30; i++)
        {
            output_30y[i] = output_30y[i-1] + L / 30;
        }
        for(int i = 0; i < 30; i++)
        {
            double f = output_30y[i] - L - R / r;
            if(f == 0)
                f = 1;
            output_30x[i] = (output_30y[i] - L) / f;
        }
    }

    void xy6(double *output_2x, double *output_2y, double k, double Dn, double Dk, double U0, double Ub)
    {
        if(Dn == Dk)
            throw exception();
        if(k == 0)
            k = 1;
        if(Ub == 0)
            Ub = 1;
        output_2x[0] = ((U0 + k * Dn) * 100) / Ub;
        output_2x[1] = ((U0 + k * Dk) * 100) / Ub;
        output_2y[0] = Dn;
        output_2y[1] = Dk;
    }

    void x4(double *input_output_x, int x_size, double d, double k)
    {
        if(k == 0)
            k = 1;
        for(int i = 0; i < x_size; i++)
        {
            input_output_x[i] = (input_output_x[i] - d) / k;
        }
    }

    void x3(double *input_df, double *output_x, int gh_size, double f0)
    {
        for(int i = 0; i < gh_size; i++)
        {
            output_x[i] = (33684210.528 / (f0 + input_df[i])) - 1402.508772;
        }
    }

    void x2(double *input_output_x, int x_size, double u1, double u2, double v1, double v2, double Ub)
    {
        double v_dif = v2 - v1;
        if (v_dif == 0)
            v_dif = 1;
        if (Ub == 0)
            Ub = 1;
        double A = (u2 - u1) / v_dif;
        double B = u2 - A * v2;
        for(int i = 0; i < x_size; i++)
        {
            input_output_x[i] = (input_output_x[i] * A + B) * 100 / Ub;
        }
    }

    void xy7(double *input_x1, double *input_y1, int gh1_size, double *input_x2, double *input_y2, int gh2_size, double *output_7poly)
    {
        Approx *gh = new Approx();
        gh->least_squares_arr(input_y1, input_x1, gh1_size, 7);
        double *temp_x1 = new double[gh2_size];
        gh->calc_polynom_func(input_y2, temp_x1, gh2_size);
        gh->least_squares_arr(temp_x1, input_x2, gh2_size, 7);
        gh->get_poly(output_7poly);
        delete[] temp_x1;
    }

    void xy8(double *input_x1, double *input_y1, double *input_x2, double *input_y2, int gh_size, double *output_7poly)
    {
        Approx *gh = new Approx();
        gh->least_squares_arr(input_y2, input_x2, gh_size, 7);
        double *temp_x2 = new double[gh_size];
        gh->calc_polynom_func(input_x1, temp_x2, gh_size);
        gh->least_squares_arr(temp_x2, input_y1, gh_size, 7);
        gh->get_poly(output_7poly);
        delete[] temp_x2;
    }
}
